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Review  of  Work  To  Date 

Our  work  falls  natiu^ly  into  three  categories:  i)  local  density  calculations  of  cluster  struc¬ 
tures  and  the  development  of  classical  potentials,  (ii)  simulations  of  thermodynamic  properties  of 
clusters  and  their  fusion,  and  (iii)  macroscopic  (continuum  theory)  modeling  of  the  sintering  pro¬ 
cess.  These  topics  will  be  addressed  in  turn.  Although  the  nominal  focus  of  the  proposal  is  on 
sintering  processes,  many  issues  important  to  sintering  are  generic  in  nature  and  thus  also  impor¬ 
tant  in  other  areas  of  cluster  science.  These  issues  have  been  approached  from  a  broad  perspec¬ 
tive,  leading  to  results  which  impact  the  whole  field  of  clusters  and  cluster-assembled  materials. 

A.  Local  Density  Calculations,  Quantum  Nfolecular  Dynamics,  and  the  Development 
of  Classical  Potentials  ' 

The  objectives  of  the  local  density  calculations  and  simulations  are  to  i)  obtain  reliable  in¬ 
formation  for  the  structure  and  energetics  of  small  and  medium  size  clusters;  (ii)  develop  meth¬ 
ods  for  first  principles  simulations  of  clusters;  (iii)  establish  contact  with  the  experimental  com¬ 
munity  by  calculating  observables  and  comparing  to  the  experimental  data  (cluster  structures  are 
not  yet  directly  measurable);  and  (iv)  develop  potentials  for  long-time  classical  simulations  of 
large  clusters. 

We  have  concentrated  on  Al,  which  is  a  "typical"  free  electron  metal,  and  Cu,  as  a 
representative  transition  metal  atom.  For  Al,  we  have  carried  out  extensive  Car-Parrinello  and 
DVM  calculations  for  clusters  containing  up  to  55  atoms.  Larger  clusters,  containing  up  to  2000 
atoms,  were  studied  quantum-mechanically  using  a  jellium  model.  For  Cu,  DVM  calculations  for 
clusters  containing  up  to  13  atoms  have  served  as  a  basis  for  improving  existing,  efficient  but 
very  approximate  “embedded  atom  method”  (EAM)  potentials.  The  results  of  the  calculations 
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were  used  to  interpret  experimental  data  obtained  by  R.  E.  Smalley  and  R.  L.  Whetten,  who  are 
also  participants  in  the  ONR  cluster  initiative. 


The  quantum  molecular  dynamics  pan  of  this  project  was  also  supported  by  the  Pittsburgh 
and  the  North  Carolina  Supercomputer  Centers,  whose  Peer  Review  Boards  have  granted  541 
and  120  hours,  respectively,  for  the  projects  described  below. 

In  our  studies  of  the  structures  of  A1  clusters  we  have  chosen  to  focus  on  large  clusters, 
since  the  motivation  for  the  present  work  is  the  smdy  of  clusters  as  building  blocks  for  the  pro¬ 
duction  of  cluster-assembled  materials.  Obviously,  the  structures  of  the  clusters  will  play  in  an 
important  role  in  the  assembly  process  and  thereby  also  influence  the  properties  of  the  final  ma¬ 
terial.  The  structural  models  most  commonly  used  for  metal  clusters  include  fragments  of  bulk 
lattices  as  well  as  icosahedral  structures.^*  ^  Icosahedral  structures  are  considered  because  they 
are  believed  to  be  the  lowest  energy  structures  for  noble  gas  clusters  up  to  very  large  cluster 
sizes.^  Since  metallic  bonding  is  non-directional  and  results  in  highly  coordinated  structures,  it 
has  been  widely  believed  that  icosahedra  are  the  ground  state  structures  for  metallic  clusters.  In 
particular,  early  theoretical  calculations  for  atoms  interacting  with  a  spherically  syiiunetric  pair 
potential  resulted  in  a  prediction  that  the  icosahedral  structures  will  be  prefnred  for  up  to  ~  2000 
atoms."*  However,  Montano  et  al*  concluded  from  the  analysis  of  EXAFS  data  that  19-55  atom 
Cu  clusters  assume  fee  structures.  On  the  other  hand,  the  authors  of  Ref.  (2)  concluded  from  the 
analysis  of  STEM  data  that  Pd  clusters  smaller  than  20  A  are  icosahedral.  These  analyses  suggest 
that  this  transition  can  occur  much  earlier  for  real  metals  and  that  its  onset  may  strongljy^pend 
on  the  kind  of  atoms  present  in  the  cluster. 


Our  calculations  for  A1  clusters  were  carried  out  using  both  the  Car-Parrine 

DVM®  methods.  The  current  formulation  of  the  Car-Parrinello  method  works  most  efficiently 

when  plane  waves  are  used  as  a  basis  set.  This  requires  the  use  of  periodic  boundary  conditions 
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and  the  clusters  were  embedded  in  a  large  simple  cubic  unit  cell  (ranging  from  30  to  40  bohr,  de¬ 
pending  on  the  size  of  the  cluster).  Due  to  the  size  of  our  clusters  we  had  to  optimize  the  pseu¬ 
dopotential  in  order  to  obtain  converged  results  with  relatively  small  plane  wave  cutoffs.  Tests 
have  shown  that  the  structural  energy  differences  calculated  here  are  well-converged  with  respect 
to  both  the  supercell  size  and  the  cutoff.'^  The  DVM  calculations  have  used  a  numerical  basis  set 
consisting  of  the  atomic  core  orbitals  and  the  3s,  3p,  and  4s  orbitals.^ 

The  smallest  clusters  for  which  both  perfect  cuboctahedral  and  icosahedral  structures  exist 
contain  13,  19,  and  55  atoms.  The  corresponding  structures  are  shown  in  Fig.  1.  Although  the 
bulk  interatomic  distance  is  a  good  first-order  guess  for  the  interatomic  distances  in  a  cluster,  the 
geometry  clearly  needs  to  be  optimized.  In  the  initial  calculations  we  have  considered  the  ener¬ 
getics  of  the  two  competing  structures,  allowing  only  breathing  mode  relaxations.  In  Fig.  2,  the 


Fig.  1  The  structures  of  the  13-,  19-,  and  55-atom  cuboctahedra  and  icosahedra. 
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binding  energy  difference  between  the  icosahedral  and  the  cuboctahedral  structures  is  shown  for 
two  plane  wave  cutoffs;  4  and  6  Ry.  It  is  evident  from  the  figure  that  the  icosahedron  is  preferred 
for  the  13-atom  cluster,  but  that  the  cuboctahedron  is  slightly  preferred  over  the  double  icosahe¬ 
dron  for  the  19-atom  cluster.  For  the  55-atom  cluster,  the  energy  difference  between  the  cubocta¬ 
hedron  and  the  Mackay  icosahedron  increases,  even  on  a  per  atom  basis.  DVM  calculations  car¬ 
ried  out  at  the  Univmity  of  Chicago  show  also  that  the  undistorted  13-atom  icosahedron  is  more 
stable  that  the  cuboctahedron  and  that  this  order  reverses  for  the  55-atom  clusters.^ 

These  initial  results  thus  predict  a  transformation  to  an  fcc-based  structure  at  ~  55  atoms.  In 
order  to  investigate  this  further  we  attempted  to  relax  all  the  structural  parameters  of  the  55-atom 
clusters  using  steepest  descent  and  conjugate  gradient  algorithms  without  imposing  symmetry 
constraints.  For  the  icosahedral  cluster,  however,  the  rate  of  convergence  of  these  methods  was 
impractically  slow.  A  slight  heating  of  the  cluster,  which  gave  the  atoms  initial  velocities,  fol¬ 
lowed  by  a  molecular  dynamics  relaxation  with  a  friction  term,  found  the  minimum  at  a  rate  of 
20-30  times  quicker  than  the  straight-forward  minimization.  This  was  due  to  the  fact  that  the 
atoms  largely  kept  their  velocities  during  the  relaxation,  thereby  traversing  quickly  through  even 


Fig.  2  The  binding  energy  difference  between  the  two  alternative  fully  symnwtric  struc¬ 
tures:  the  icosahedron  and  the  cuboctahedron.  Breathing  mode  relaxations  are  included. 
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Table  I.  Structural  energy  differences  relative  to  the  radially-relaxed  Mackay  icosahedron  for  the 
55-atom  A1  clustn^  [eV/clusterJ.  The  plane  wave  cutoff  was  4  Ry  (6  Ry). 


radially  relaxed  j 

"weakly-annealed " 

■  ITf? 

cuboctahedron 

-2.1  (-1.9) 

-5.3 

the  flatter  parts  of  the  potential  energy  surface. 

The  "weak  annealing"  procedure  resulted  in  the  reversal  of  the  energies  of  the  two  struc¬ 
tures,  with  the  icosahedron  now  being  preferred  by  over  1  eV  over  the  cuboctahedron  (see  Table 
I).  One  should  note,  however,  that  since  the  energies  of  the  two  structures  remained  within  ~  1 
eV,  while  the  relaxation  energies  ranged  from  3  to  6  eV,  a  search  for  a  global  minimum  was  nec¬ 
essary.  Due  to  the  size  of  these  systems,  however,  this  task  was  far  from  trivial.  After  considering 
several  alternatives  we  adapted  the  efficient  "constant  thermodynamic  speed  annealing"  method 
to  molecular  dynamics.^  Our  tests  for  a  number  of  clusters  have  shown  that  this  method  is  sub¬ 
stantially  more  reliable  than  competing  methods  in  finding  the  global  minimum  from  a  variety  of 
starting  configurations.  In  what  is  probably  to  date  the  largest  quantum  system  studied  by  simu¬ 
lated  annealing,  we  have  succeeded  in  transforming  an  initially  cuboctahedral  55-atom  cluster  to 
an  icosahedron.  Annealing  of  the  icosahedron,  on  the  other  hand,  does  not  change  its  structure. 
These  results  confirm  that  the  "weakly-annealed"  icosahedron  is  indeed  the  lowest  energy  struc¬ 
ture  for  the  55-atom  cluster.  The  energy  difference  between  the  icosahedral  and  the  cuboctahe¬ 
dral  structures  is  still  unexpectedly  small  and  appears  dominated  by  surface  reconstruction  (see 
below).  The  transformation  to  fcc-based  structures  is  thus  still  expected  to  occur  very  early  for  A1 
clusters. 

Since  the  structure  of  small  clusters  cannot  be  measured  directly  at  present  and  only  indi¬ 
rect  information  is  available,  we  set  out  to  calculate  Ionization  Potentials  (IPs)  and  Electron 
Affinities  (EAs)  for  these  clusters.  This  requires  the  evaluation  of  the  total  energy  for  charged 
states.  However,  the  Car-Parrinello  calculations  proceed  in  supercell  geometry  and  the  energy  of 
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a  periodic  array  of  charges  diverges.  Therefore,  we  had  to  develop  special  procedures  for  these 
calculations.  We  proceeded  by  carrying  out  CP  supercell  calculations  for  the  charged  clusters 
while  adding  a  neutralizing  background  charge  to  the  system.  The  neutralizing  charge  was  spread 
uniformly  throughout  the  entire  cell  and  it  had  a  very  small  effect  on  the  total  charge  density.  The 
CP  charge  density  without  the  neutralizing  background  was  used  for  the  construction  of  the 
Coulomb  potential  for  an  isolated  cluster,  subject  to  zero  boundary  conditions  at  infinity.  The 
whole  procedure  was  tested  by  increasing  the  size  of  the  cell,  which  lowers  the  density  of  the 
background  charge  and  improves  the  treatment  of  the  boundary  conditions  at  infinity.  The  results 
were  both  stable  and  accurate. 

The  advantage  of  the  present  procedure  for  calculating  IPs  and  EAs  is  that  the  plane  wave 
basis  set  is  used  for  the  calculations  of  charged  states.  It  has  more  variational  freedom  than 
atomic  basis  sets  which  are  optimized  for  ground  state  calculations.  The  disadvantage  is  that  the 
charge  density  used  for  the  evaluation  of  the  total  energy  of  charged  clusters  is  not  strictly  self- 
consistent.  However,  since  the  total  energy  expression  is  variational,  this  error  is  small,  which  is 
also  demonstrated  by  the  insensitivity  of  the  results  to  the  size  of  the  CP  cell  (once  it  is  suffi¬ 
ciently  large). 

The  remaining  errors  in  the  calculations  are  due  to  the  finite  size  of  the  basis  set,  the  ne¬ 
glect  of  the  spin  polarization  effects,  and  to  the  LDA.  In  order  to  further  test  the  accuracy  of  our 
procedures  as  well  as  of  the  LDA,  we  calculated  the  ionization  potential  of  the  A1  atom  using  the 
method  described  above  (with  the  same  plane  wave  cutoffs)  and  by  numerically  integrating  the 
all-electron  LDA  equations  with  and  without  spin-polarization  effects.  The  LDA  all-electron  re¬ 
sult  of  6.05  eV  (including  spin)  is  in  excellent  agreement  with  the  experimental  IP  of  5.98  eV. 
Without  spin-polarization,  the  all-electron  result  of  5.85  eV  is  in  very  good  agreement  with  the 
CP  value  of  5.77  eV.  The  difference  between  the  spin-polarized  and  the  spin-averaged  results  can 
also  be  used  as  an  estimate  of  the  upper  bound  of  the  neglect  of  spin-polarization  in  our  results 
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for  clusters.  We  believe  that  these  are  the  Tvst  high  precision  calculations  of  IPs  and  EAs  for 
clusters. 


The  calculated  and  measured  IPs  and  EAs  are  summarized  in  Table  n.  In  comparing  the 
calculated  results  to  experiment,  one  should  bear  in  mind  that  LDA  is  far  more  accurate  in  its 
predictions  of  IPs  than  EAs  for  small  clusters.  For  large  clusters  the  accuracy  should  be  similar, 
since  in  the  infinite  cluster  limit  both  the  IP  and  the  EA  become  equal  to  the  woiic  function, 
which  is  well  predicted  by  the  LDA.  Table  11  also  shows  that  the  calculated  results  are  higher 
than  experiment  by  0.4  -  1.0  eV.  The  experimental  accuracy  is  better  than  0.15  eV  for  the  abso¬ 
lute  values,  while  the  differences  between  IPs  or  EAs  for  different  cluster  sizes  are  determined 
far  more  accurately.  The  discrepancies  between  theory  and  experiment  are  partially  due  to  the 
use  of  the  LDA  approximation,  in  particular  for  the  case  of  EAs,  and  partially  due  to  the  discrep¬ 
ancies  between  the  calculated  and  measured  structures.  When  the  clusters  undergo  simulated  an- 


Table  11.  A  comparison  of  the  calculated  and  measured  IPs  and  EAs  for  several  cluster  structures 
[eV].  The  1-atom  case  is  treated  as  a  test  of  the  numerical  accuracy  of  the  Car-Parrinello  calcula¬ 
tions.  See  text. 


Exp  IP 

ExpEA 

1-atom 

numerical  LD 
numerical  LD  with  spin 
Car-Parrinello 

5.85 

6.05 

5.77 

3:5s 

3-atom  cluster 
triangle 

6.66 

6.4i  -  ^.45  (a) 

2.05 

1.5  (b),  1.5^  (c) 

13-atom  cluster 
icosahedron 
cuboctahedron 
annealed  icosahedron 

7.15  6.0 

6.62 

6.56 

mu 

3.81  2.9 

3.22 

3.04 

i.(J  (b),  i.66  (c) 

19-atom  cluster 
icosahedron 
cuboctahedron 

6.01 

5.81 

3.05 

2.81 

55-atom  cluster 
icosahedron 
cuboctahedron 
annealed  icosahedron 

5.10  5.15 

5.50  5.0 

5.35 

BNIIIIIII 

2.87  3.2 
3.29  2.9 
3.19 

(a)  from  Ref  (10);  (b)  from  Ref  (11);  (c)  from  Ref  (12);  (d)  from  Ref  (13). 


nealing,  the  agreement  between  theory  and  experiment  improves  substantially.  For  the  larger 
clusters,  structural  differences  for  the  unannealed  structures  are  indeed  the  largest  source  of  error 
(cf.  Table  H). 

An  analysis  of  the  "weakly-annealed"  structures  (Fig.  3)  has  shown  substantial  differences 
as  well  as  similarities.  In  particular,  although  the  geometric  structures  are  quite  different,  the 
structure  factors  and  radial  distribution  functions  are  very  similar.  These  results  cast  substantial 
doubt  on  the  value  of  line  shape  analysis  of  scattering  data  as  a  tool  for  the  determination  of  the 
structure  of  small  particles.  Although  not  very  apparent  from  the  figure,  the  changes  in  the  core 
of  the  clusters  upon  annealing  are  small  and  most  of  the  atomic  rearrangement  occurs  on  the  sur¬ 
face,  as  expected. 

The  analysis  of  the  A1  clusters  can  be  carried  much  further  by  considering  the  structure 
factors,  angular  distributions,  and  pair  correlation  factors  for  the  ideal,  weakly  annealed,  and 
fully  annealed  structures.  An  interesting  aspect  of  annealing  studies  which  can  be  afforded  on  to¬ 
day's  supercomputers  is  that  each  QMD  run  results  in  a  significantly  different  structure.  An 
analysis  based  on  the  structure  factors  and  distribution  functions  shows^^  that  the  local  topology 
of  each  structure  is  almost  the  same,  i.e.,  the  nearest  neighbor  and  the  most  probable  angle  peaks 
occur  at  exactly  the  same  positions  for  the  three  annealed  clusters,  and  there  are  strong  similari¬ 
ties  in  the  positions  of  second  neighbor  peaks.  It  follows  that  the  energetics  of  Alss  is  nuunly  de¬ 
termined  by  the  first  and  possibly  second  nearest  neighbors.  This  result  provides  a  justification 
for  the  success  of  the  effective  medium  and  embedded  atom  methods, in  which  the  screened 
interatomic  interactions  in  metals  are  approximated  by  short  range  effective  potentials. 

The  multitude  of  structures  for  the  55-atom  cluster  is  possible  because  nx>st  of  the  atoms 
reside  on  the  surface  and  are  therefore  free  from  the  constraints  of  bulk  packing.  Due  to  the  size 
of  this  cluster,  slight  changes  in  the  positions  of  distant  neighbors  can  lead  to  completely  in- 
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equivalent  structures.  However,  a  careful  analysis  of  the  annealed  structures  on  graphics  work¬ 
stations  capable  of  manipulating  three-dimensional  images  shows  that  the  positions  of  the  inner 
atoms  diffor  substantially  between  these  structures  so  that  the  clusters  cannot  simply  be  modeled 
as  a  common  core  and  differently  arranged  surfaces. 

We  have  also  carried  out  a  number  of  Car-Parrinello  calculations  in  order  to  establish  a  li¬ 
brary  of  structure-energy  relationships  for  the  fitting  of  classical  potentials.  When  these  results 
were  compared  to  those  obtained  using  a  bulk  EAM  potential,  we  found  that  although  the 
amount  of  radial  relaxation  was  predicted  very  well  for  a  given  geometry,  the  energy  differences 
between  different  structures  were  predicted  less  well,  and  sometimes  resulted  in  a  wrong  order¬ 
ing  of  the  structures  (Fig.  4).  We  have  considered  several  schemes  for  the  development  of  a  po¬ 
tential  specifically  fitted  to  clusters.  The  method  which  worked  best  to  date,  however,  involved 
just  fitting  the  values  of  the  embedding  function  and  its  derivative  to  the  results  of  quantum-me¬ 
chanical  calculations  for  the  equilibrium  structures  of  clusters  consisting  of  3-12  atoms.  This 
procedure  resulted  in  a  better  determination  of  the  embedding  function  in  the  region  of  small 
densities,  without  sacrificing  the  bulk  region.  The  improvement  in  the  results  for  the  13-  and  19- 
atom  clusters  in  both  the  icosahedral  and  the  cuboctahedral  geometries  was  substantial.  However, 
the  improved  embedding  function  still  predicted  the  ideal  icosahedron  to  be  lower  in  energy  than 


Fig.  3.  The  structures  of  the  weakly-annealed  cuboctahedron  and  icosahedron. 
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the  cuboctahedron.  Including  the  SS-atom  clusters  in  the  fit  resulted  in  a  correct  ordering  of  the 
structures  (sec  Fig.  4). 

Several  points  should  be  made  about  the  EAM  results  to  date:  (i)  Although  small  clusters 
were  used  in  the  fit,  only  modest  improvements  are  seen  in  the  3-8  atom  cluster  range.  This  is 
most  likely  due  to  the  fact  that  the  quantum  effects  are  the  strongest  in  this  range,  due  to  the  large 
spacing  between  the  energy  levels  in  small  clusters.  When  this  spacing  decreases,  the  EAM  de¬ 
scription  should  get  better.  For  cluster  sizes  greater  than  12,  the  improvement  is  substantial;  (ii) 
T!ie  fitting  of  the  55-atom  results  has  resulted  in  a  variation  of  the  curvature  of  the  embedding 
function,  although  the  function  .tself  has  remained  monotonic;  and  (iii)  The  distorted  geometries 
resulting  from  energy  minimization  in  the  quantum-mechanical  calculations  lead  to  an  increase 
in  the  EAM  total  energy  even  when  the  improved  embedding  function  is  used.  Due  to  the  com¬ 
plexity  of  the  potential  energy  surfaces  uncovered  by  the  Car-Parrinello  calculations,  we  believe 
that  further  improvements  in  this  potential  are  needed. 

For  clusters  larger  than  ~  1(X)  inequivalent  atoms  even  the  fastest  methods  for  quantum- 
mechanical  calculations  for  clusters  become  too  expensive.  However,  spheroidal  clusters 
containing  thousands  of  atoms  can  still  be  studied  utilizing  a  jellium  approximation,  which  is 


Clustar  size 


Fig.  4.  The  difference  between  the  quantum-mechanical  and  EAM  results,  using 
the  bulk  A1  EAM  potential  and  the  cluster-derived  (improved)  EAM  potential. 
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known  to  work  very  well  for  bulk  Al.  Indeed,  the  work  described  above  has  shown  that  the 
structural  changes  occur  mainly  at  the  surface,  and  in  very  large  clusters  inner  atoms  out-number 
surface  atoms  by  a  substantial  margin.  Our  jellium  calculations  have  been  carried  out  for  clusters 
containing  up  to  2000  atoms,  and  found  that  they  exhibit  groupings  of  electronic  energy  levels 
into  what  have  been  called  supershells.  Such  groupings  have  been  sf  m  in  the  experiments  of 
our  collaborator,  R.  L.  Whetten.  In  contrast  to  the  levels  of  the  large  clusters,  we  have  found  by 
comparison  to  results  of  our  DVM  calculations®  that  the  jellium  model  does  not  represent  well 
the  electronic  level  structure  of  smaller  clusters,  e.g.  AI55.  The  size,  boundary  effects  and 
effective  crystal  field  together  cause  the  level  pattern  to  deviate  strongly  from  the  pattern  implied 
by  the  jellium  model.  The  convergence  to  jellium-like  behavior  may  become  a  useful  index  of 
the  extent  of  metal-like  behavior  in  small  clusters. 

We  have  also  initiated  the  construction  of  an  embedded  atom^5  Qj  potential  for  clusters, 
utilizing  previous  EAM  results  as  well  as  our  DVM  calculations  for  several  structures  of  Cu6, 
Cu9,  and  CU13.  Our  strategy  here  is  to  adjust  the  parameters  of  the  simple  potentials  to  fit  results 
of  local  density  calculations  in  the  vicinity  of  the  most  important  points  on  the  multidimensional 
potential  surface,  namely  minima  and  saddle  points.  This  procedure  is  now  feasible  because  lo¬ 
cating  saddle  points  is  now  almost  as  easy  as  locating  minima.  To  date,l®  the  adjustments  of 
the  embedded  atom  potential  have  only  been  made  in  the  regions  of  the  secondary  minima,  and 
not  yet  in  the  regions  of  the  saddle  points.  However,  the  feasibility  of  the  method  has  been  estab¬ 
lished  and  its  full  implementation  is  one  of  the  tasks  for  our  next  phase. 

Although  the  Car-Parrinello  method  is  quite  general,  its  most  efficient  implementation  re¬ 
lies  on  the  use  of  a  plane  wave  basis  set.  Periodic  boundary  conditions  are  necessarily  imposed, 
and  calculations  involving  first-row  or  transition  metal  atoms  become  quite  expensive,  some¬ 
times  prohibitively  so.  In  order  to  reduce  the  computational  effort  for  such  clusters,  we  have  de¬ 
veloped  a  new  real  space  method  based  on  multigrid  techniques  which  can  be  used  with  or  with- 
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out  periodic  boundary  conditions.  It  can  also  handle  first-row  and  transition  metal  atoms  without 
a  substantial  increase  in  computational  time.  In  this  method,  the  Kohn-Sham  and  Poisson's  equa¬ 
tions  are  solved  on  several  grids  at  a  time.  The  final  solution  is  obtained  on  a  sufficiently  fine 
grid  so  that  it  rq)resents  well  the  pseudo-wavefunctions  and  the  pseudopotentials  in  the  problem. 
However,  the  convergence  is  accelerated  by  the  use  of  several  auxiliary  grids,  each  with  the 
spacing  increased  by  a  factor  of  two  over  the  previous  grid.  Denoting  the  grid  spacing  by  Ax,  one 
can  show  that  iterations  on  any  given  grid  reduce  quickly  those  Fourier  components  of  the  error 
with  wavelengths  ~  Ax  but  are  ineffective  when  dealing  with  substantially  longer  wavelengths. 
The  role  of  the  extra  grids  is  thus  to  converge  quickly  and  systematically  all  wavelength  compo¬ 
nents  of  the  error.  One  important  advantage  of  the  algorithm  is  that  the  grids  need  not  be  uni¬ 
formly  spaced  so  that  the  grids  can  be  enhanced  locally  near,  e.g.,  transition  metal  atoms. 
Substantial  programming  and  some  algorithm  development  efforts  are  still  needed,  but  we  al¬ 
ready  obtained  test  results  for  H,  H2,  Li,  and  LiH.  Some  of  these  results  have  been  described  in 
invited  talks  at  the  Faraday's  bicentennial  ^9  and  the  cluster  conference  in  Richmond.20 

A  major  new  effort  initiated  under  this  grant  and  continued  in  a  separate  grant  is  a  study  of 
buckyballs  and  related  structures.  Under  this  grant,  we  computed  the  equilibrium  atomic  struc¬ 
ture  of  solid  Qjo  studied  its  dynamics  during  heating  to  high  temperatures.^!  The  calcula¬ 
tions  determined  the  precise  atomic  positions  of  the  carbon  atoms  (which  were  much  later  con¬ 
firmed  by  neutron  scattering),  showed  that  C^o  start  to  rotate  at  relatively  low  temperatures  even 
when  in  a  solid  form  (which  explained  NMR  data),  and  that  they  are  stable  even  at  very  high 
temperatures.  These  results  received  a  very  wide  exposure  and  were  quoted  in  several  scientific 
and  popular  periodicals.  The  electron  distribution  figures  and  snapshots  from  the  QMD  simula¬ 
tions  have  appeared  on  the  covers  of  Science  News,  Science,  Scientiflc  &  Engineering  Indicators 
1992,  and  in  several  books  and  foreign  periodicals. 


B.  Classical  Molecular  Dynamics  Simulations 

Our  simulations  have  focused  on  the  issues  associated  with  meidng  and  surface  melting  of 
clusters,  which  have  important  consequences  for  cluster  fusion  and  sintering;  on  modeling  of 
cluster  fusion;  and  on  the  development  of  an  analytical  quantum-statistical  theory  of  surface 
melting.  The  last  project  has  been  to  a  large  extent  guided  by  the  results  of  the  simulations, 
which  is  the  reason  for  its  inclusion  into  this  section.  As  in  the  case  of  our  quantum  calculations, 
the  issues  have  been  approached  from  broad  perspective,  resulting  in  generic  results  with  wide 
applicability. 

The  simulations  of  single  clusters  have  shown  that  clusters  of  fewer  than  about  45  atoms 
show  no  behavior  that  could  be  called  surface  melting.  Instead,  some  display  simple  melting,  and 
others  show  softening,  manifested  by  passage  among  a  limited  set  of  catchment  basins  on  the 
cluster’s  multidimensional  surface,  and  then,  at  higher  energies,  a  transition  to  liquid-like  behav¬ 
ior.  This  means  that  "liquid-phase  sintering",  the  process  by  which  the  necks  of  sintering  parti¬ 
cles  are  filled  by  liquid-like,  mobile  material,  is  a  concept  applicable  to  clusters  containing  two 
or  more  shells  of  atoms.  However  for  clusters  of  about  50  or  more  atoms,  surface  melting  seems 
to  be  the  general  rule,  occurring  over  a  range  of  temperatures  well  below  the  bulk  melting  tem¬ 
perature.  Figure  5  illustrates  the  diffusion  coefficients  in  melted  surfaces,  as  functions  of  the 
mean  cluster  temperature,  for  Lennard-Jones  clusters  of  four  sizes. 

The  criteria  for  surface  melting  are  the  same  as  those  for  bulk  melting  of  clusto^,  but  can 
be  applied  only  by  designating  the  particles  in  the  cluster's  surface  to  be  distinguishable  from 
those  in  the  core  of  the  cluster.  Making  this  distinction  is  straightforward  to  incorporate  into  sim¬ 
ulations  by  classical  mechanical  molecular  dynamics  (MD),  the  procedure  we  have  used  for  most 
of  the  simulations  for  this  work.  These  criteria  include:  the  "Lindemann  criterion"  of  large-ampli¬ 
tude  motions  of  nearest  neighbors,  meaning  10%  or  more  of  the  equilibrium  interparticle  dis- 
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tance;  diffusion,  typically  evaluated  from  the  mean  square  displacement  of  the  atoms  from  their 
initial  locations,  as  functions  of  time,  (which  is  a  necessary  but  not  sufficient  condition  for  liq¬ 
uid-like  behavior);  the  occurrence  of  vibrational  modes  of  very  low  frequency  ("soft  modes”, 
also  a  necessary  but  not  sufficient  condition);  passage  of  the  system  among  many  potential  wells 
at  a  rate  not  very  different  from  the  frequency  of  vibrations  within  one  well,  which  is  probably 
both  necessary  and  sufficient  if  the  set  of  wells  visited  includes  those  of  permutational  isomors  of 
the  structure  of  lowest  energy;  and  passage  among  all  the  permutational  isomers  of  the  structure 
of  lowest  energy,  which  is  a  sufficient  but  not  a  necessary  condition  for  liquid-like  behavior. 
There  are  others,  such  as  pair  or  radial  distribution  functions  and  angular  distribution  functions 
(which  give  necessary  but  not  sufficient  conditions  for  liquid-like  behavior  because  amorphous 
solids  look  like  liquids  according  to  these  criteria),  which  we  also  sometimes  use. 
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Fig.  6.  Examples  of  snapshots  of  two  (Ar)i47  clusters  during  the  process  of  sintering:  a) 

after  1000  steps  of  10’*5  sec;  b)  after  10,000  steps;  c)  after  40,000  steps;  d)  after  100  000 
steps. 

Simulations  of  Lennard-Jones  clusters  ranging  in  size  from  45  to  147  atoms  have  provided 
diffusion  coefficients  for  these  clusters  as  functions  both  of  temperature  and  of  cluster  size. 
Examples  are  shown  in  Fig.  5.  While  the  diffusion  coefficients  for  clusters  of  fixed  size  are 
smooth  functions  of  temperature,  they  are  not  smooth  functions  of  cluster  size,  at  least  in  the 
lower  end  of  our  size  range.  It  will  be  necessary,  in  order  to  make  the  general  theory  readily  ap¬ 
plicable  to  clusters  of  all  important  sizes  for  sintering,  to  carry  the  simulations  of  surface  melting 
to  clusters  large  enough  that  the  size  dependence  of  the  diffusion  coefficient  is  smooth.  This  will 
be  a  direct  continuation  of  work  already  done,  and  will  not  be  particularly  demanding,  although  a 
substantial  amount  of  computer  time  will  be  needed. 
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The  simulation  of  sintering  itself  has  been  carried  out  for  pairs  of  Leonard- Jones  clusters  of 
147  atoms.  Such  clusters  exhibit  a  combination  of  liquid-phase  sintering  and  reconstruction 
which  allows  an  approximate  separation  of  the  process  into  an  early  stage  of  neck  formation  and 
a  second  stage  of  cementing  and  reconstruction,  during  which  the  two  clusters  still  retain  some 
visible  identity.  Copies  of  photographs  of  some  of  the  computer  animations  of  these  simulations 
are  shown  in  Fig.  6. 

Simulations  of  surface  melting  on  single  clusters  of  147  argon  (Lennard- Jones)  atoms  re¬ 
vealed  unforeseen  characteristics  of  that  process.  Individual  frames  showing  the  cluster  at  arbi¬ 
trary  instants  along  its  MD  trajectory  fh  our  prior  expectations;  the  core  of  55  atoms  retains  its 
identity  as  a  polyhedron  with  only  small-amplitude  distortions  but  the  atoms  of  the  surface  seem 
to  be  amorphous,  with  no  obvious  polyhedral  structure.  However,  movies  of  the  same  simula¬ 
tions  show  that  the  real  situation  is  not  just  the  fonnation  of  a  liquid  layer  around  a  solid.  In  the 
energy  range  within  which  the  diagnostics  say  that  the  surface  is  molten  and  the  core  is  solid,  the 
situation  is  rich  and  complex.  Most  of  the  core  atoms,  in  any  brief  time  interval,  execute  large- 
amplitude,  anharmonic  oscillations  about  equilibrium  positions  on  the  polyhedron  of  the  solid 
cluster,  while  one  or  a  few  atoms  play  the  role  of  "floaters",  moving  rather  easily  around  the  sur¬ 
face  of  the  cluster.  Over  time,  the  "floaters"  exchange  roles  with  the  other  atoms  of  the  surface 
layer,  so  that  permutations  and  diffusion  are  achieved.  The  set  of  surface  atoms  meets  all  the  cri¬ 
teria  of  liquid-like  behavior,  but  by  microscopic  means  quite  different  from  what  was  tradition¬ 
ally  expected  of  liquids.  On  this  basis  it  is  important  to  reexamine  the  thermodynamic  analysis  of 
surface  melting  of  clusters,  which  had  previously  implied  that  liquid  surfaces  of  clusters  are  sta¬ 
ble  only  as  fluctuations  of  thermally  equilibrated  solid  clusters.^ 

To  give  a  basis  to  the  simulations  of  surface  melting,  we  developed  an  analytic,  micro¬ 
scopic  quantum-statistical  theory  of  surface  melting.  This  is  parallel,  in  some  ways,  to  the  tiieory 
we  developed  for  homogeneous  melting  of  clusters,^^’^^  but  has  a  different  form  precisely  be* 
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cause  of  what  we  have  learned  about  "floaters"  and  "stay-at-homes"  from  the  simulation  movies. 
A  combination  of  simulation  and  analytic,  statistical-mechanical  theoiy  led  to  a  new,  quantifiable 
model  for  surface  melting  of  clusters  and  to  a  set  of  diffusion  coefficients  for  clusters  of  various 
sizes  at  various  temperatures.  A  full  report  of  this  work  in  now  in  press.^^  The  conclusions  of 
this  work  are  that  surface  melting  of  clusters  of  atoms  bound  by  central  forces  is  a  real  phe¬ 
nomenon  for  clusters  of  more  than  about  50  atoms,  that  the  mechanism  of  this  melting  is  primar¬ 
ily  the  promotion  of  a  few  surface  atoms  from  the  outer  shell  to  a  "floating"  state  outside  the 
normal  outermost  shell,  leaving  a  vacancy  there,  that  most  of  the  diffusive  motion  iz  carried  by 
the  "floaters"  which  exchange  occasionally  with  atoms  in  the  surface  layer,  that  the  surface¬ 
melting  transition  has  only  a  single  stable  form  at  any  temperature  and  pressure  and  is  in  this 
sense  like  a  second-order  phase  transition  rather  than  a  first-order  transition  (as  homogeneous 
melting  of  clusters  is),  and  that  because  the  process  involves  formation  of  floaters,  it  has  a  small 
latent  heat  and  is  therefore  not  a  true  second-order  phase  transition. 

A  second  achievement  of  this  group  was  the  development,  in  collaboration  with  Professor 
Robert  Whetten  and  his  student  Xiuling  Li,  of  a  full  phase  diagram  for  a  cluster.  (Professor 
Whetten's  work  is  being  carried  out  under  the  same  ONR  Initiative  as  our  own.)  This  is  the  first 
time  that  full  diagrams  of  the  curves  of  equal  mean  chemical  potential,  in  terms  of  pressure  and 
volume,  pressure  and  temperature,  and  volume  and  temperature,  as  well  as  heat  capacity  and 
thermal  expansion  coefficients,  have  been  developed  for  a  realistic  cluster  of  atoms.  This  work  is 
now  in  press^^  and  is  also  being  presented  at  the  Spring  Meeting  of  the  Materials  Research 
Society  in  San  Francisco,  1992. 

In  order  to  address  nonmetallic  clusters  relevant  to  real  applications,  particularly  to  sinter¬ 
ing  of  refractory  powders,  we  have  investigated  the  behavior  of  salt  clusters,  particularly  (KCl)n. 
These  were  chosen  because  they  can  be  modeled  reliably  and  are  similar  in  many  ways  to  metal 
oxide  clusters,  which  are  of  obvious  general  interest.  Ousters  of  NaO  have  been  studied  for  a 
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variety  of  reasons,2**29  but  not  for  the  purposes  we  followed:  surface  melting,  the  possibility  of 
glass  formation  and,  ultimately,  sintering.  The  results  to  date  have  not  yet  given  a  clear  answer 
about  surface  melting,  largely  because  we  have  not  yet  simulated  clusters  large  enough  to  be  ex¬ 
pected  to  show  that  phenomenon.  The  results  regarding  glass  formation  are  striking:  it  appears 
that  clusters  as  small  as  (KCl)i6  are  capable  of  being  quenched  to  glasses.30  The  stability  of  these 
glasses  is  not  yet  known,  and  will  be  one  of  the  most  important  questions  to  be  addressed  in  the 
future.  Bulk  ionic  glasses  have  been  recognized  for  some  time.3132 

C.  Continuum  Modeling 

The  continuum  theory  part  has  achieved  several  of  its  initial  goals.  We  have  developed  an 
analytic  theory  of  the  shapes  of  sintering  bodies  in  first-stage  sintering,33, 34, 35, 36  a  quantitative 
theory  of  the  unpinning  and  elimination  of  grain  boundaries  in  sintering^?,  and  a  quantitative 
theory  of  third  stage  sintering38,  in  the  arena  of  simulations,  an  intensive  database  of  random 
close  packed  structures  of  distributions  of  sphere  sizes  has  been  amassed  and  used  to  report  and 
model  the  statistics  of  contact  distributions^^.  These  distributions  were  used  to  simulate  early 
stage  sintering  to  study  the  relative  effects  of  grain  boundary  versus  surface  transport^O*^!.  The 
computations  associated  with  this  part  of  the  effort  amounted  to  over  500  Cray  hours.  This  time 
was  provided  by  San  Diego  Supercomputer  Center  via  a  block  grant  to  SDSU. 

Since  the  driving  force  for  sintering  is  the  reduction  in  surface  and  grain  boundary  areas 
during  particle  fusion,  our  study  of  macroscopic  models  began  with  the  observation  that  the  static 
geometric  structure  of  the  sinter,  representing  the  solid  skeleton  and  the  mobile  pool,  may  be 
characterized  as  a  new  version  of  the  problem  of  minimal  surfaces:  Plateau's  problem.  The  new 
version  follows  simply  from  the  exotic  nature  of  the  "liquid"  pool  which  covers  the  solid  skele¬ 
ton.  The  exoticness  of  this  liquid  derives  from  the  near  indistinguishability  of  the  solid  and  liquid 
forms  making  the  surface  tension  between  liquid  and  solid  zero  and  the  surface  tensions  of  \apar 
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against  solid  or  vapor  against  liquid  equal.  Funhermore,  microscopic  pore  sizes  make  gravita¬ 
tional  effects  negligible.  This  leaves  us  with  a  much  simplified  special  case  of  the  classical  the¬ 
ory  of  capillarity.  This  case  was  considered  unphysical  and  left  by  the  wayside  by  previous  gen¬ 
erations  of  mathematicians  and  physicists.  It  is  however,  the  physically  appropriate  model  for 
sintering  processes.  This  discovery  has  already  lead  to  some  interesting  developments  in  capil¬ 
larity  and  in  understanding  the  geometry  of  microporous  materials  and  is  certain  to  lead  to  many 
more.  It  turns  out  that  the  special  case  lends  itself  to  a  different  mathematical  formulation^^ 
which  is  much  more  amenable  to  analytic  attack.  This  allowed  us  to  construct  nearly  analytic 
solutions  for  the  early  stages  of  sintering.  ^<33  It  also  allowed  us  to  prove  such  theorems  as:  The 
mean  curvature  on  a  solid  portion  of  the  skeleton  surface  is  everywhere  greater  than  the  mean 
curvature  of  the  liquid  portion.  Perhaps  most  importantly,  it  promises  a  universality  of  structure 
among  pores  found  in  sinters.  Our  formulation  of  the  generalized  Plateau  problem  represents  an 
important  spinoff  application  that  this  project  has  produced  for  the  theory  of  minimal  surfaces. 
The  generalized  problem  leads  naturally  to  several  other  interesting  new  problems  in  the  theory 
of  minimal  surfaces^2^ 

Our  analysis  beyond  the  static  structure  has  leaned  heavily  on  the  assumption  that  the  rate 
limiting  step  for  the  sintering  process  is  the  creation  of  a  mobile  pool  of  material.  For  conve¬ 
nience,  this  pool  will  be  called  "liquid"  below  although  the  only  requirement  in  the  formalism  is 
that  the  mobility  of  this  pool  be  much  greater  than  the  diffusion  rate  of  the  atoms  in  the  core.  The 
latter  assumption  has  been  shown  to  hold  by  microscopic  simulations  (see  preceding  subsection). 
As  a  consequence  of  this  assumption,  the  redistribution  of  the  "liquid"  will  lead  to  an  equilibrium 
configuration  that  minimizes  the  surface  area,  which  is  proportional  to  the  free  enei^  of  the  sys¬ 
tem.  The  local  mean  curvature  represents  a  generalized  force  responsible  for  the  motion  of  the 
mobile  pool.  It  follows  that  the  mean  curvature  of  any  optimal  surface  is  constant  wherever  it  is 
covered  by  the  mobile  pool,  since  such  surfaces  represent  equilil»ium  configurations. 
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The  generation  of  optimal  surfaces  did  not  turn  out  to  be  straightforward  however.  We  be¬ 
gan  with  the  case  of  two  identical  solid  spherical  particles  in  point  contact,  wetted  by  a  certain 
volume  of  liquid.  The  assumption  of  spherical  particles  is  mild  since  any  spheroidal  shaped  par¬ 
ticles  become  spherical  once  a  mobile  pool  of  atoms  is  created.  The  situation  for  identical  two 
spheres  was  described  by  a  first  order  differential  equation,  which  turned  out  to  be  singular  in  an 
important  region  of  the  parameter  space.  Schemes  were  found  for  stabilizing  the  singularity  and 
the  equation  was  finally  solved  by  numerical  integration.  The  integrand  contained  two  Lagrange 
parameters  whose  values  had  to  be  matched  by  shooting  methods.  An  overlap  parameter  was 
added  which  (if  negative)  could  also  describe  the  case  of  separated  spheres.  The  generated  solu¬ 
tions  could  simulate  the  sintering  of  two  spheres  and  of  a  close  packed  configuration  of  spheres 
with  any  desired  ratio  of  surface  to  grain  boundary  melting.  Furthermore,  these  simulations 
evolved  the  solid  skeleton  as  well  as  the  mobile  liquid  pool  and  thus  went  far  beyond  the  original 
solutions  of  the  static  case.  The  results  for  two  identical  spheres  were  published  in  Ref.  (34) 
where  it  was  seen  that  in  some  cases  significant  differences  existed  between  the  shapes  calcu¬ 
lated  by  minimizing  the  free  energy  and  the  shapes  assumed  by  previous  workers. 

The  simulations  above  used  our  analytic  and  numerical  work  on  the  problem  of  unequal 
spheres.  New  methods  were  developed  including  an  explicit  algorithm  for  the  construction  of  the 
minimal  surfaces  needed  to  solve  the  sintered  shape  problem  in  any  axially  symmetric  casc.3536 
We  noted  there  that  the  problem  is  fundamentally  combinatorial  in  nature  with  many  cases  de¬ 
termined  by  the  geometry  of  the  solid  to  be  wetted.  We  showed  that  these  cases  could  be  conve¬ 
niently  handled  by  noting  the  volumes  of  certain  "separatrix"  solutions,  which  restrict  the  range 
of  choices  for  the  boundary  conditions.  The  identification  of  these  separatrix  solutions  is  impor¬ 
tant  for  maintaining  stability  and  convergence  properties  of  the  algorithms. 

The  global  problem  of  minimizing  the  riee  energy  by  optimally  allocating  the  liquid  among 
the  puddles  is  again  a  difficult  combinatorial  optimization  problem.  Our  discovery,  that  it  is  pos- 
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sible  to  map  this  problem  into  a  dynamic  programming  problem^  is  an  important  innovation 
both  for  the  simulations  in  this  study  and  the  understanding  of  microporous  geometry. 

For  unequal  spheres,  the  first  observation  was  that  the  point  of  minimal  constriction  be¬ 
tween  two  unequal  size  spheres  moves  in  response  to  the  force  to  maintain  the  smallest  possible 
surface  area  as  liquid  is  added.  Once  the  constriction  between  the  spheres  disappears,  the  grain 
boundary  becomes  unpinned  and  anneals  out.  We  have  carried  out  essentially  exact  simulations 
of  grain  boundary  annealing.  The  first  surprise  was  that  the  unpinning  always  occurs  when  the 
liquid  neck  between  the  two  spheres  is  exactly  tangent  to  the  top  of  the  smaller  sphere.  The  sec¬ 
ond  surprise  was  that  the  thinnest  point  in  the  neck  does  not  necessarily  move  directly  toward 
this  point;  its  initial  movement  is  toward  the  center  of  the  larger  sphere,  provided  that  grain 
boundary  transport  is  fast  compared  to  surface  transpon,  i.e.  the  primary  source  of  the  mobile 
pool  is  the  grain  boundary .37 

Generating  the  solutions  for  unequal  size  spheres  in  terms  of  the  physically  important  vari¬ 
ables  such  as  liquid  volume  and  sphere-sphere  overlap  (rather  than  the  Lagrange  parameters) 
turned  out  to  be  also  a  numerically  demanding  problem.  In  order  to  simplify  most  of  the  time 
evolution  runs  we  decided  to  create  a  database  of  solutions  suitable  for  splining.  This  data  base 
was  assembled  and  will  prove  invaluable  for  the  continuing  work: 

Database  1:  Numerical  solutions  of  the  shapes  (and  all  geometric  parameters)  of  the  opti¬ 
mal  surfaces  required  to  describe  the  wetting  of  two  spheres  of  radii  rj  and  T2  with  separation  h 
using  an  anoount  Vl  of  liquid  volume.36<3S  Using  this  information  we  can  give  a  complete  de¬ 
scription  of  the  evolution  of  two  solid  spherical  particles  for  any  ratio  of  surface/grain  boundary 
transport. 33  Furthermore,  the  comparison  between  neighboring  optimal  surfaces  for  infinitesimal 
changes  in  the  anangement  of  the  solid  skeleton  gives  the  forces  acting  on  the  solid  cores  of  the 
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particles.  The  solutions  thus  provide  nearly  the  complete  information  necessary  fen-  simulating 
the  time  evolution  of  the  structure  of  the  sinter. 


Our  next  concern  was  to  start  with  a  reasonably  accurate  model  of  the  sinter.  Accordingly, 
we  generated  a  second  database: 

Database  2:  Configurations  of  random  packed  spheres  and  disks  with  various  distributions 
of  radii.^^  This  data  base  has  provided  stadstics  for  modeling  the  distribution  of  sphere-sphere 
contacts.  It  also  provides  the  initial  configurations  for  simulating  sintering  processes. 

The  simplest  model  for  time-evolution  of  the  sinter  uses  the  random  packings  in  database  2 
to  provide  the  number  and  type  of  sphere-sphere  contacts  and  evolves  these  independently  em¬ 
ploying  the  infonnation  stored  in  database  1.  This  is  a  direct  application  of  the  methods  in  Ref. 
(34)  for  a  distribution  of  sphere  sizes  and  necks.  We  have  also  explored  versions  of  this  model 
which  allow  partial  interaction  between  neighboring  necks  by  allocating  mobile  volume  among 
the  necks  in  an  optimal  fashion.^3 

While  such  models  give  qualitatively  useful  information  -  sec  Amar  et  al.^  for  the  case  of 
identical  spheres  and  Schdn  et  al.^>  for  the  case  of  a  random  close  packed  array  of  spheres  ~ 
they  neglect  the  change  with  time  of  the  three  dimensional  structure  of  the  solid  material  beyond 
the  independent  pairs  approximation.  They  also  neglect  the  possible  coalescence  of  neighboring 
necks. 

The  worit  described  above  deals  with  early  stage  sintering;  we  have  also  had  good  success 
in  late  stage  sintering.  Following  standard  usage,  we  distinguish  three  different  stages  of  the  sin¬ 
tering  process.^  During  first  stage  sintering,  the  nnobile  material  forms  Inidges  between  neigh¬ 
boring  sintering  particles,  and  fills  in  narrow  necks.  During  second  stage  sintering,  the  channels 
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that  connect  various  regions  of  the  sinter  close,  creating  weakly  connected  pores  in  the  process. 
These  pores  will  become  approximately  spherical  at  the  conclusion  of  this  stage.  Finally,  in  the 
third  stage,  the  pores  exchange  material  via  connecting  grain  boundaries,  leading  to  growth  of  Me 
bigger  pores  at  the  expense  of  the  smaller  ones;  this  is  the  so-called  coarsening  process. 

For  third  stage  sintering  the  situation  is  simple  enough  for  an  analytic  attack.  We  have 
modeled  the  third  stage  as  a  coarsening  process  of  the  Lifshitz-Slyozov  type^^  that  proceeds  via 
grain  boundary  diffusion.^^  The  slow  time  scale  is  still  the  promotion  of  a  particle  to  the  "mobile 
pool,"  but  in  this  case  the  mobile  pool  is  just  the  material  in  the  grain  boundary  connecting  dif¬ 
ferent  spherical  pores.  Since  the  distances  between  the  isolated  pores  are  large,  each  pore  effec¬ 
tively  equilibrates  with  the  reservoir  of  material  in  the  grain  boundaries,  which  in  turn  serves  as 
the  medium  for  the  interaction  among  the  pores  leading  to  a  growth  of  the  larger  pores  at  the  ex¬ 
pense  of  the  'mailer  ones.  Thus  far  we  are  dealing  with  the  standard  Lifshitz-Slyozov  model.^5 
The  distinction  comes  in  observing  that  the  grain  boundary  represents  only  a  two  dimensional 
medium  coupling  our  three  dimensional  pores.  We  solved  the  Lifshitz-Slyozov  model^5  for  the 
case  of  an  m  dimensional  medium  surrounding  n  dimensional  pores. The  specific  example  of 
m=2,  n=3  was  treated  in  detail  to  deduce  the  time  evolution  of  the  size  distribution  of  pores  con¬ 
nected  by  grain  boundaries,  as  well  of  the  average  pore  size  and  the  number  of  pores  per  unit 


Fig.  7.  Three  competing  equilibrium  structures  for  filling  a  pore  between  four  identical  disks. 
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volume. 


The  problem  of  second  stage  sintering  is  much  more  difficult  to  simulate.  The  reason  is  al¬ 
ready  evident  in  two  dimensions,  where  much  more  is  accessible  analytically.  As  an  illustration, 
consider  a  pore  between  four  identical  disks.  This  simple  case  already  reveals  an  unexpected 
complexity  of  structure^^  with  three  competing  local  minima:  symmetric,  asymmetric,  and 
bridged  (see  Rg.  7).  Each  of  these  types  represents  the  energetically  favored  configuration  over 
some  range  of  the  area  parameter. 

Even  ignoring  the  problems  due  to  local  minima  such  as  these,  the  interaction  between 
necks  presents  several  new  difficulties,  which  are  common  to  both  the  two-  and  three-dimen¬ 
sional  case.  They  are  associated  with  keeping  track  of  which  necks  have  coalesced  and  how  to 
classify,  lookup,  and  calculate  the  associated  "many-neck"  puddles.  Since  we  are  able  to  investi¬ 
gate  these  issues  more  easily  for  two  dimensional  sintering,  we  plan  to  begin  exploring  second 
stage  sintering  by  implementing  the  simulation  in  two  dimensions.  The  program  and  associated 
data  structures  handle  two  or  three  dimensional  problems  with  equal  ease  so  the  transition  to 
three  dimensional  data  is  anticipated  to  be  relatively  easy.  The  two  dimensional  case  will  develop 
the  necessary  intuition  for  modeling  the  more  involved  three-dimensional  case.  It  is  also  ex¬ 
pected  to  reveal  the  appropriate  universality  classes  of  pore  structures  in  two  dimensions.  The 
transition  to  three  dimensions  can  be  achieved  using  the  Surface  Evolver  package  developed  by 
the  Minnesota  Supercomputer  Center.  Preliminary  work  implementing  this  software  has  shown 
that  it  is  capable  of  the  task.  The  calculation  for  three  spheres  is  presently  being  generated. 
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